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We present an efficient strategy for controlling a vast range of non-integrable quantum many- 
body one-dimensional systems that can be merged with state-of-the-art tensor network simulation 
methods like the Density Matrix Renormalization Group. To demonstrate its potential, we employ 
it to solve a major issue in current optical-lattice physics with ultra-cold atoms: we show how to 
reduce by about two orders of magnitudes the time needed to bring a superfluid gas into a Mott 
insulator state, while suppressing defects by more than one order of magnitude as compared to 
current experiments [J]. Finally, we show that the optimal pulse is robust against atom number 
fluctuations. 



Classical control theory has played a major role in the 
development of present-day technologies [2| . Likewise, re- 
cently developed quantum optimal control methods [3|— 15|] 
can be applied to emerging quantum technologies, e.g. 
quantum information processing - until now, at the level 
of a few qubits However, such methods encounter 

severe limits when applied to many-body quantum sys- 
tems: due to the complexity of simulating the latter, 
existing quantum control algorithms (requiring many it- 
erations to converge) usually fail to yield a desired fi- 
nal state within an acceptable computational time. A 
paradigmatic application of control of many body quan- 
tum system is the control of the dynamics of a quantum 
phase transition. The process of crossing a phase transi- 
tion in an optimal way has been studied for decades for 
classical systems. Only recently it has been recast in the 
quantum domain, attracting a lot of attention (see e.g. Q 
and references therein) since, for instance, it has implica- 
tions for adiabatic quantum computation and quantum 
annealing 10]. A transition between different phases is 
usually performed by "slowly" (adiabatically) sweeping 
an external control parameter across the critical point, 
allowing for a transformation from the initial to the fi- 
nal system ground state with sufficiently high probabil- 
ity. However, at the critical point in the thermodynami- 
cal limit a perfect adiabatic process is forbidden in finite 
time Thus, the resulting final state (for finite-time 
transformations) is characterized by some residual exci- 
tation energy, corresponding to the formation of topolog- 
ical defects within finite-size domains. The Kibble-Zurek 
theory has been shown to yield good estimates of the 
density of defects or of the residual energy 0, [ljj . The 
importance of these estimates in the quantum domain 
is underscored by the fact that, apart from very specific 
cases where analytical solutions are available, theoreti- 
cal investigations must rely on heavy numerical simula- 
tions due to the exponential growth of the Hilbert space 
with the system size and to the diverging entanglement 
at the critical point [l3j. Nevertheless, it is possible to 



perform one-dimensional simulations of the dynamics by 
means of tensor-network-based techniques such as the 
time-dependent Density Matrix Renormalization Group 
(t-DMRG) 0. 

The basic underlying idea of classical control is to pick 
a specific path in parameter space to perform a specific 
task. This is formally a cost functional extremization 
that depends on the state of the system and is attained by 
varying some external control parameters. In a quantum- 
mechanical context, a big advantage is that the goal can 
be reached via interference of many different paths in pa- 
rameter space, rather than just one. In few-body quan- 
tum systems, it has been shown that optimal control finds 
optimal paths in the parameter space that result in con- 
structive interference of the system's classical trajectories 
toward a given goal [Bj]. Indeed, present optimal control 
strategies demonstrated an impressive control of quan- 
tum systems, ranging from optimization of NMR pulses 
8] to atomic [H| and superconducting qubits [f| , as well 
as the crossing of a quantum phase transition (OPT) in 
the analytically solvable quantum Ising model 16]. How- 
ever, despite their effectiveness, they cannot be efficiently 
applied to systems that require tensor network methods 
for their simulation. 

The present letter marks a step further — it provides 
for the first time a means to control the evolution of a 
non-integrable many-body quantum system, resulting in 
the optimization of a given figure of merit. This is done 
by introducing a strategy to integrate optimal control 
with t-DMRG simulations of the many-body quantum 
dynamics. 

Tensor-network-simulations- Tensor network methods 
are based on the assumption that it is possible to describe 
approximately a wide class of states with a simple ten- 
sor structure. In particular the DMRG describes ground 
states static properties of one dimensional systems by 
means of a Matrix Product State (MPS) [T3]- The main 
characteristic of an MPS is that the resources needed 
to describe a given system depend only polynomially on 
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FIG. 1: A) An initial guess pulse c°(t) is used as a starting 
point. B) The function J- (a) for the case a = {a\,a,2\ and 
the initial polytope (white triangle) are defined and moved 
"downhill" (darker red triangles) until convergence is reached. 
C) The final point is recast as the optimal pulse c{i). 



the system size N, due to the introduction of an an- 
cillary dimension m that determines the precision of the 
approximation. Since an exact description requires expo- 
nentially increasing resources with the number of compo- 
nents N, the tensor network approach results in an expo- 
nential gain in resources. Given a system Hamiltonian, 
the best possible approximated description of the system 
ground state -within the MPS at fixed m- is determined 
by means of an efficient energy minimization. With some 
slight modification, discretizing the time T = n steps Af 
and performing a Trotter expansion, the algorithm can 
be adapted to follow a state time evolution, the so-called 
t-DMRG [Hj]. The t-DMRG is a very powerful numeri- 
cal method for efficiently numerically simulating the time 
evolution of one-dimensional many body quantum sys- 
tems. The class of states and of time evolutions that can 
be efficiently described with a small error is determined 
by the presence of entanglement between the different 
system components [l3[. Here, we will use the t-DMRG 
for the simulation of cold atoms in time dependent optical 
lattices, which we feed into the Chopped RAndom Basis 
(CRAB) optimization algorithm as described below. 

The CRAB method- The general scenario of an op- 
timal control problem can be stated as follows: given 
a system described by a Hamiltonian H depending on 
some control parameters Cj(t) with j — 1, . . . ,Nc, the 
goal is to find the c/s time dependence (pulse shape) 
that extremizes a given figure of merit T , for instance 
the final system energy, state fidelity, or entanglement. 
We then start with an initial pulse guess c°(t) and look 
for the best correction that has a simple expression in a 
given functional basis. As an explicative example, here 
we focus on the case where the correction is of the form 
Cj(t) = c9 (t) • fj(t), and the functions fj (t) can be simply 
expressed in a truncated Fourier space, depending on the 
expansion coefficients <Xj = ah [k = 1, . . . , Mj). In par- 
ticular, in the following, we start from an initial ansatz, 
e.g. an exponential or linear ramp, and we introduce a 
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FIG. 2: Scheme of the Mott-superfluid transition in the ho- 
mogeneous system for average occupation number (n) = 1: 
increasing the lattice depth V (black line) the atoms super- 
fluid wave functions (upper) localize in the wells (lower). If 
the transition is not adiabatic -or optimized- defects appear 
(here represented by a hole and a doubly occupied site). 



correction of the form 



/(*) 



Af 



l + ^Ak sin(i/fct) + B k cos{v k t) 



(1) 



Here, k — 1, . . . , M, v k = 27rfc(l + rk)/T are "random- 
ized" Fourier harmonics, T is the total time evolution, 
rfc G [0 : 1] are random numbers with a flat distribution, 
and Af is a normalization constant to keep the initial 
and final control pulse values fixed. The optimization 
problem is then reformulated as the extremization of a 
multivariable function J-({Ak), {B^}, {vk)}, which can 
be numerically approached with a suitable method, e.g., 
steepest descent or conjugate gradient [l8j]. When using 
CRAB together with t-DMRG, computing the gradient 
of J- is extremely resource consuming, if not impossible. 
Thus we resort to a direct search method like the Nelder- 
Mead or Simplex methods [11]. They are based on the 
construction of a polytope defined by some initial set of 
points in the space of parameters that "rolls down the 
hill" following predefined rules until reaching a (possibly 
local) minimum (see Fig. [T]). Due to the fact that direct 
search methods are based on many independent evalua- 
tions of the function to be minimized, they can be effi- 
ciently implemented together with t-DMRG simulations 
(and possibly performed in parallel). We stress that the 
functional dependency of the correction presented here 
(Eq. (JT])) is one possible approach: different strategies 
might be explored. Indeed, making a given choice con- 
fines the search of the optimal driving field in a subspace 
of the whole space of functions and the results might de- 
pend on this choice. On the other hand, this approach 
simplifies the optimization problem that would be other- 
wise computationally unfeasible when t-DMRG simula- 
tions are needed. As shown below, the described choice 
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FIG. 3: Initial guess (dashed black) and optimal ramp (solid 
red) V(t) for the Bose Hubbard model in the presence of the 
trap with N = 30 sites, total time evolution T ~ 3 ms. Inset: 
Populations (m) (empty black) and fluctuations (Anf) (full 
red) at time t = T for the exponential initial guess (circles) 
and optimal ramp (squares) for N = 10. 



allows to perforin a successful optimization. 

The optical-lattice system- Very recently, the experi- 
mental and theoretical analysis of the dynamics of cold 
atoms in optical lattices has experienced a fast develop- 
ment, after the experimental demonstration of coherent 
manipulation of ultra-cold atoms in the seminal work of 
Ref. [l9( , where a Bose-Einstein condensate is first loaded 
into a single trap, and then a periodic lattice potential is 
slowly ramped up, inducing a quantum phase transition 
to a Mott insulator. This is the enabling step for a wide 
range of experiments, from transport or spectroscopy to 
quantum information processing [20 ] . In most of these 
applications, it is essential to achieve the lowest possible 
number of defects in the final state, that is, to reach ex- 
actly a final state with fixed number of atoms per site, 
e.g. unit filling. Up to now, this has been pursued by 
limiting the process speed - the superfluid-Mott insula- 
tor transition has been performed in about a hundred 
ms, with a density of defects typically of the order of 
10% |H|. 

Cold atoms in an optical lattice can be described by the 



Bose Hubbard model defined by the Hamiltonian [20ll26| 



H= 



3 



N , 



-j(b]b j+1 +h.c.)+n(j-- 



(2) 

The first term on the r.h.s. of Eq.((5]) describes the tun- 
neling of bosons between neighbouring sites with rate J; 
fi is the curvature of the external trapping potential, and 
rij = bjbj is the density operator with bosonic creation 

(annihilation) operators &j (bj) at site j = 1, . . . , N. The 
last term is the on-site contact interaction with energy 
U. The system parameters U and J can be expressed as 
a function of the optical lattice depth V (we set h — 1 
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FIG. 4: Residual defect density p for N = 10,20,30,40, 
T ~ 3ms, p c = 0.001 for the homogeneous system (green 
squares) and in the presence of the trap (grey circles). The 
red region highlights the typical p for different initial ramp 
shapes (see text). Inset: Final p computed applying the pulse 
optimized for system size N = 20 to different system sizes 
A7V = —4, ... ,4 (at constant filling). The results are almost 
independent from the truncation error for m > 50. 



from now on) 20]. As sketched in Fig. [2j the system 
undergoes a quantum phase transition from a supcrfluid 
phase to a Mott insulator as a function of the ratio J/U. 
In a homogeneous one-dimensional system, the QPT is 
expected to occur at J c /U ~ 0.083, where (upon decreas- 
ing the ratio J/U) the ground state wave function dras- 
tically changes from a Fermi- Thomas distribution with 
high fluctuations in the number of particles per site to a 
simple product of local Fock states with no fluctuations 
in the number of atoms per site [20I ] . In the presence of an 
external trapping potential on top of the optical lattice, 
the phase diagram is more complex: the two phases coex- 
ists in different trap regions and typical "cake" structures 
are formed 12211. 



Results- Following previous numerical studies [23J that 
modeled the experiment [l[ , and supported by strong ev- 
idence of agreement between numerical simulations and 
experimental results 0, H^] , we studied both the ideal 
ho mog eneous system (O = 0) and the experimental setup 
of [24J where the trapping potential is present. We ap- 
plied the CRAB optimization to the preparation of a 
Mott insulator with ultra-cold atoms in an optical lat- 
tice, that is, we optimized the ratio J/U(t) that drives 
the superfluid-Mott insulator transition. The resulting 
optimal ramp shape drives the system into a final Mott 
insulator state with a density of defects below half a per- 
cent in a total time of the order of a few milliseconds, 
amounting to a drastic improvement in the process time 
and in the quality of the final state - by about two orders 
of magnitude and by more than one, respectively. 
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We consider a starting value of the lattice depth 
V(0) = 2E r corresponding to J/U(0) ~ 0.52, since 
the description of the experimental system by (Eq. [2]) 
breaks down for U(0) < 2E r (26|. However, the ini- 
tial lattice switching on (V = — > 2E r ) can be per- 
formed very quickly without exciting the system (few 
milliseconds at most) [27j . We optimize the ramp to 
obtain the minimal residual energy per site AE/N — 
(E(T)—Eg)/N (where Eg is the exact final ground state 
energy). In all simulations performed we set the total 
time T = 50h/U ~ 3.01ms and the final lattice depth 
V(T)/E r = 22 ~ 2.4 ■ 10~ 3 J/U, deep inside the Mott in- 
sulator phase. Unless explicitly stated, we set the average 
occupation to one (J2i( n i) = ^0- ln a ^ DMRG simula- 
tions, we exploited the conservation of number of parti- 
cles and used m = 20, . . . , 100, At = 10~ 2 -r- 10~ 3 . We 
computed the final density of defects p = jjJ2i \( n i) ~M : 
when it reached a given threshold p c — 10 -3 , the opti- 
mization was halted. In Fig. [3] we report a typical result 
of the optimization process: the initial guess and final op- 
timal ramp for the system in the presence of the confining 
trap are shown for the parameter values corresponding 
to the experiment [l|, for a system size N = 30. As it 
can be clearly seen, the pulse is modulated with respect 
to the initial exponential guess and no high frequencies 
are present, reflecting the constraint introduced by the 
CRAB optimization. In the inset we display the final 
occupation numbers and the corresponding fluctuations, 
for the initial exponential ramp and the optimal pulse 
in the case N — 10. The figure clearly demonstrates 
the convergence to a Mott insulator in the latter case as 
fluctuations are drastically reduced and the occupation 
is exactly one for every site. 

Finally, in Fig. 0] we plot the optimized density of de- 
fects p as a function of the system size (up to N = 40) 
for the homogeneous and for the trapped system, demon- 
strating an improvement with respect to the initial guess 
by one (two) orders of magnitudes. Indeed, the expo- 
nential guess - like other guesses: linear, random, and a 
pulse optimized for a smaller system (TV = 8 sites) - gave 
residual density of defects of the order of 10% (red region 
in Fig. 01). To check the experimental feasibility of our 
findings, we studied the stability of the optimal evolu- 
tions under different sources of error and experimental 
uncertainties, like atom number fluctuations. The in- 
set of Fig. [4] shows the final density of defects when an 
optimal pulse computed for a given system size, is ap- 
plied to a different system size (keeping the average fill- 
ing constant). As it can be seen, the optimization works 
also for system size fluctuations of up to 20%: the fi- 
nal density of defects is of the same order. This robust- 
ness is crucial as the experimental realization of these 
systems is performed in parallel on many different one- 
dimensional tubes with different numbers of atoms (20j . 
We also checked the cases of different filling and of pulse 
distortion obtaining similar results (data not shown). 



Outlook- In conclusion, we would like to mention that 
the CRAB optimization strategy introduced here can in 
principle be applied also to open quantum many-body 
systems, e.g. by means of recently introduced numer- 
ical techniques 28j . Perhaps an even more stimulat- 
ing perspective would be that of implementing it with 
a quantum system in place of the t-DMRG classical sim- 
ulator, i.e. performing a CRAB based closed-loop op- 
timization [29l | . The optimization might be performed 
during the experimental repetitions of the measurement 
processes, thus adding a small overhead to the experi- 
mental complexity. This would extend the applicability 
of the CRAB method to the optimization of quantum 
phenomena that are completely out of reach for simula- 
tion on classical computers, and represent a major design 
tool for future quantum technologies. 

We thank F. Dalfovo, J. Denschlag, R. Fazio, C. Fort, 
M. Greiner, C. Koch, C. Menotti, G. Pupillo and M. Rizzi 
for discussions; the PwP project (www.dmrg.it), the 
DFG (SFB/TRR 21), the bwGRiD, the EC (grant 
247687, AQUTE, and PICC) for support. 



[1 
[2 

[3 

[1 
[5 

[6 

[7: 

[8 
[9 

[10 

[11 
[12 

[13 

[11 
[15 
[16 
[17 

[is: 

[19 
[20 

[21 
[22 

[23 
[24 



T. Stoferle, et. al, Phys. Rev. Lett. 92, 130403 (2004). 
J.T. Betts "Practical Methods for Optimal Control using 
Nonlinear Programming", SIAM, Philadelphia, (2001). 

D. D'Alessandro "Introduction to Quantum Control and 
Dynamics", Chapman & Hall/CRC (2007). 

C. Brif, R. Chakrabarti, and H. Rabitz arXive:0912.5121. 
V. F. Krotov, "Global Methods in Optimal Control The- 
ory" (Marcel Dekker, New York, 1996). 

S. Montangero, T. Calarco, and R. Fazio, Phys. Rev. 
Lett. 99, 170501 (2007). 

P. Rebentrost, et. al, Phys. Rev. Lett. 102, 090401 
(2009). 

N. Khaneja, et. al, J. Magn. Reson. 172, 296 (2005). 
W.H. Zurek, U. Dorner, P. Zoller Phys. Rev. Lett. 95 
105701 (2005). 

E. Farhi, et. al, Science 292 472 (2001). 

M. Born and V.A. Fock, Zeit. fur Physik 51 165 (1928). 
W.H. Zurek, Phys. Rep. 276, 177 (1996). 
G. Vidal, Phys. Rev. Lett. 91, 147902 (2003). 
Schollwock, Rev. Mod. Phys. 77, 259 (2005). 
Calarco, et. al, Phys. Rev. A 70 , 012306 (2004). 
Caneva, et. al, llarXiv: 1011.66341 

Ostlund and S. Rommer, Phys. Rev. Lett. 75, 3537 
(1995); J. I. Cirac and F. Verstraete, Journ. Of Phys. A 
42, 504004 (2009). 

W.H. Press, et. al, "Numerical Recipes", Cambridge 

University Press, New York 2007. 

M. Greiner, et. al, Nature 415, 39 (2002). 

I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 

80, 885 (2008). 

D. Jaksch, P. Zoller, Annals of Physics, 315, 52 (2005). 
G.G. Batrouni, et. al, Phys. Rev. Lett. 89, 117203 
(2002). 

C. Kollath, et. al, Phys. Rev. Lett. 97, 050402 (2006). 
CD. Fertig et. al, Phys. Rev. Lett. 94 120403 (2005). 



■5 



[25] I. Danshita and C.W. Clark, Phys. Rev. Lett. 102, 

030407 (2009). 
[26] D. Jaksch, et. at, Phys. Rev. Lett. 81, 3108 (1998). 
[27] M.Greiner, private communication. 



[28] M. Zwolak, G. Vidal, Phys. Rev. Lett. 93, 207205 (2004). 
[29] C. Brif, R. Chakrabarti, and H. Rabitz, New Journ. of 
Phys., 12 075008 (2010). 



